!> author: 左志华
!> date: 2022-07-04
!> version: alpha
!>
!> Smoothed kernel function <br>
!> 光滑核函数
module sph_smoothed_kernel_function
    use sph_kinds, only: rk
    use seakeeping_constants, only: Pi
    real(rk), parameter :: r23 = 2.0_rk/3.0_rk  !! 2/3
    !> 光滑核函数
    type smoothed_kernel_function
        real(rk) :: selfden  !! 粒子的自身密度加权系数
        real(rk) :: hsml  !! 光滑长度
        integer :: scale = 2  !! 光滑长度缩尺比例
        real(rk) :: scale_hsml  !! 缩尺后的光滑长度
        character(:), allocatable :: type  !! 光滑函数类型
        procedure(kernel_fcn), pointer :: kernel => cubic_spline_kernel  !! 光滑核函数句柄
        real(rk), private :: factor  !! 核函数常量
    contains
        procedure :: setup_kernel
    end type smoothed_kernel_function
    type(smoothed_kernel_function) :: skf
    abstract interface
        !> Kernel function <br>
        !> 光滑核函数
        pure subroutine kernel_fcn(skf, r, dx, w, dwdx)
            import :: rk, smoothed_kernel_function
            class(smoothed_kernel_function), intent(in) :: skf  !! Smoothed kernel function <br>
                                                                !! 光滑核函数
            real(rk), intent(in) :: r                   !! Euclidean distance between two points <br>
                                                        !! 欧氏距离
            real(rk), intent(in) :: dx(2)               !! dimensional distance between two points <br>
                                                        !! 对应坐标轴距离
            real(rk), intent(out) :: w                  !! kernel value <br>
                                                        !! 核函数
            real(rk), intent(out) :: dwdx(2)            !! derivative of kernel value <br>
                                                        !! 核函数导数
        end subroutine kernel_fcn
    end interface
contains
    !> 建立光滑核函数句柄
    subroutine setup_kernel(this, hsml)
        class(smoothed_kernel_function), intent(inout) :: this
        real(rk), intent(in) :: hsml
        real(rk) :: hv(2)
        this%hsml = hsml
        select case (skf%type)
        case ('cubic_spline')
            this%scale = 2
            this%kernel => cubic_spline_kernel
            this%factor = 15.0_rk/(7.0_rk*Pi*this%hsml*this%hsml)
        case ('quintic_spline')
            this%scale = 3
            this%kernel => quintic_spline_kernel
            this%factor = 7.0_rk/(478.0_rk*Pi*this%hsml*this%hsml)
        case default
            this%scale = 2
            this%kernel => cubic_spline_kernel
            this%factor = 15.0_rk/(7.0_rk*Pi*this%hsml*this%hsml)
            print '(a)', 'Warning: Unknown kernel type, use cubic spline kernel by default'
        end select
        call this%kernel(0.0_rk, [0.0_rk, 0.0_rk], this%selfden, hv)
        this%scale_hsml = this%hsml*this%scale
    end subroutine setup_kernel
    !> Smoothed kernel function: Cubic spline <br>
    !> 三次样条曲线 (monaghan 1985) <br>
    !> 计算光滑函数 \( W_{ij} \) 及其导数 \( \frac{dW}{dx_{ij}} \) 的子例程
    pure subroutine cubic_spline_kernel(skf, r, dx, w, dwdx)
        class(smoothed_kernel_function), intent(in) :: skf
        real(rk), intent(in) :: r           !! Euclidean distance between two points <br>
                                            !! 欧氏距离
        real(rk), intent(in) :: dx(2)       !! Dimensional distance between two points <br>
                                            !! 对应坐标轴距离
        real(rk), intent(out) :: w          !! kernel value <br>
                                            !! 核函数值
        real(rk), intent(out) :: dwdx(2)    !! Derivative of kernel value <br>
                                            !! 核函数导数
        real(rk) :: q
        q = r/skf%hsml
        if (q <= 1.0_rk) then
            w = skf%factor*(r23 - q*q + q**3/2)
            dwdx = skf%factor*(-2.0_rk + 1.5_rk*q)*dx/(skf%hsml*skf%hsml)     ! 计算机原则，先乘后除
        else if (q > 1.0_rk) then
            w = skf%factor*(2.0_rk - q)**3/6.0_rk
            dwdx = -skf%factor*(2.0_rk - q)**2*dx/(2.0_rk*skf%hsml*r)
        end if
    end subroutine cubic_spline_kernel
    !> Smoothed kernel function: Quintic spline <br>
    !> 四次样条曲线 (Liu 2010)
    pure subroutine quintic_spline_kernel(skf, r, dx, w, dwdx)
        class(smoothed_kernel_function), intent(in) :: skf
        real(rk), intent(in) :: r           !! Euclidean distance between two points <br>
                                            !! 欧氏距离
        real(rk), intent(in) :: dx(2)       !! Dimensional distance between two points <br>
                                            !! 对应坐标轴距离
        real(rk), intent(out) :: w          !! kernel value <br>
                                            !! 核函数值
        real(rk), intent(out) :: dwdx(2)    !! Derivative of kernel value <br>
                                            !! 核函数导数
        real(rk) :: q
        q = r/skf%hsml
        if (q <= 1.0_rk) then
            w = skf%factor*((3 - q)**5 - 6*(2 - q)**5 + 15*(1 - q)**5)
            dwdx = skf%factor*(-120.0_rk + 120.0_rk*q - 50.0_rk*q**2)*dx/skf%hsml**2
        else if (q > 1.0_rk .and. q <= 2.0_rk) then
            w = skf%factor*((3 - q)**5 - 6*(2 - q)**5)
            dwdx = skf%factor*(-5.0_rk*(3 - q)**4 + 30.0_rk*(2 - q)**4)*dx/(skf%hsml*r)
        else if (q > 2.0_rk .and. q <= 3.0_rk) then
            w = skf%factor*((3 - q)**5)
            dwdx = skf%factor*(-5.0_rk*(3 - q)**4)*dx/(skf%hsml*r)
        else
            w = 0.0_rk
            dwdx = 0.0_rk
        end if
    end subroutine quintic_spline_kernel
end module sph_smoothed_kernel_function
